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We show that the formation of a magnon condensate in thin ferromagnetic films can be ex¬ 
plained within the framework of a classical stochastic non-Markovian Landau-Lifshitz-Gilbert equa¬ 
tion where the properties of the random magnetic field and the dissipation are determined by the 
underlying phonon dynamics. We have numerically solved this equation for a tangentially magne¬ 
tized yttrium-iron garnet film in the presence of a parallel parametric pumping field. We obtain a 
complete description of all stages of the nonequilibrium time evolution of the magnon gas which is in 
excellent agreement with experiments. Our calculation proves that the experimentally observed con¬ 
densation of magnons in yttrium-iron garnet at room temperature is a purely classical phenomenon 
which should be called Rayleigh-Jeans rather than Bose-Einstein condensation. 
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In the past decade the nonequilibrium dynamics of 
parametrically pumped magnons in thin yttrium-iron 
garnet (YIG) films has been investigated by many experi¬ 
mental studies m- Very rich physics was found, includ¬ 
ing the overpopulation of the lowest energy state, which 
was interpreted as Bose-Einstein condensation (BEC) of 
magnons at room temperature and finite momentum. Us¬ 
ing the technique of Brillouin light scattering it is even 
possible to measure the magnon distribution with mo¬ 
mentum and time resolution m ■ This allows an obser¬ 
vation of the parametric resonance and of the subsequent 
thermalization leading to the formation of the conden¬ 
sate in detail mm- Unfortunately, a complete theoreti¬ 
cal understanding of this phenomenon is still lacking and 
there is no theory that can simultaneously describe all 
stages of the experiment. While the so-called S-theory 
mifTHl is able to describe the parametric resonance used 
to populate certain magnon states, it does not properly 
take magnon-magnon scattering into account and there¬ 
fore cannot describe the cascade of relaxation processes 
leading to the formation of a magnon condensate. On the 
other hand, theories that focus on the condensate usu¬ 
ally do not take the pumping dynamics into account and 
start with some given quasiequilibrium state which can 
be identified with the ground state of some effective quan¬ 
tum mechanical Hamiltonian mm- Phenomenological 
approaches of the Ginzburg-Landau type also have been 
used to study the condensation dynamics ,223- Finally, 
theories dealing with the relaxation processes and kinet¬ 
ics of excited magnons did not include the possibility of 
magnon condensation 

Since BEC is a manifestation of quantum mechan¬ 
ics, it seems at the first sight reasonable that quantized 
magnons obeying Bose statistics are essential to obtain 
a satisfactory theoretical description of magnon conden¬ 
sation in YIG. However, since the experiments are per¬ 
formed at room temperature, which is large compared 
with the relevant magnon energies, the equilibrium distri¬ 


bution of the magnons is the Rayleigh-Jeans rather than 
the Bose-Einstein distribution. This suggests that the 
experimentally observed magnon condensation at room 
temperature m is a purely classical phenomenon which 
should be called Rayleigh-Jeans condensation rather than 
BEC. The concept of Rayleigh-Jeans condensation has 
been discussed before by Kirton and Keeling in the con¬ 
text of condensation of photons m ■ In fact, it is well 
known that classical waves subject to some randomness 
and mode-coupling can undergo a kinetic condensation 
transition analogous to BEC of quantum systems (55U30] . 
Recently this phenomenon has been observed by imaging 
classical light dynamics in a photorefractive crystal [30] . 
In this Letter we show that the condensation of magnons 
in YIG HHS] is another physical realization of kinetic con¬ 
densation of classical waves. 

In order to calculate the nonequilibrium time evolu¬ 
tion of the pumped magnon gas in YIG, we model the 
spin dynamics in YIG by means of a stochastic Landau- 
Lifshitz-Gilbert m equation (LLG) with non-Markovian 
damping due to the coupling to the thermal phonon bath. 
The classical description is justified because the relevant 
magnon energies are much smaller than room tempera¬ 
ture where the experiments are performed. Moreover, the 
saturation magnetization of YIG is rather large so that 
the spins can be approximated by three-component clas¬ 
sical vectors Si(t) of length S ss 14 which are localized at 
the sites Ri of a cubic lattice [ 35 ] 133 ] ■ The coupling of the 
spins to the thermal phonon bath gives rise to random¬ 
ness and dissipation. Often these terms are taken into 
account phenomenologically in the LLG. Fortunately, for 
YIG the microscopic form of the dominant spin-phonon 
coupling is well known [M[ [35] , so that we can derive 
the LLG microscopically by eliminating the phonon co¬ 
ordinates from the equations of motion of the coupled 
spin-phonon system, as described in Refs. [36] 37] ■ In this 
way we arrive at the following stochastic LLG with non- 
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Markovian damping, 


Si(t) 


Si(t ) x 


Si(t) x 


3 

[ dt'y2<G ij (t,t')s j (t'). 
Jo j 


(1) 


Here the time-dependent external magnetic field H ( t ) is 
measured in units of energy and Ky is a matrix in the 
spin components with elements K = 5 a ^Jij + , 

where the exchange couplings Jy connect nearest neigh¬ 
bor sites with strength J, and DfJ is the dipolar tensor 
[33133133 • The effect of the phonon bath on the spin 
dynamics manifests itself via the induced magnetic field 
hi(t) and the non-Markovian dissipation kernal G ij(t, t’), 
which is a matrix in the spin components. A microscopic 
derivation of these quantities starting from the coupled 
spin-phonon Hamiltonian for YIG is given in the Supple¬ 
mental Material E3- It turns out that the dissipation 
kernel G ij(t,t') can be expressed in terms of the proper¬ 
ties of the phonon bath as follows: 



= E £ e ik ( R '- R i) 

[ns kX 

n \n .cos[u kx (t - t')] 
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where N is the number of lattice sites, uj kx is the disper¬ 
sion of an acoustic phonon with polarization vector e kXl 
and A labels one longitudinal and two transverse phonon 
polarizations. The effective ionic mass in the unit cell 
is denoted by M, and we have introduced the notation 
k a p = k a ep + kpe ai where e a are Cartesian unit vectors 
in direction a = x,y,z. Because for YIG the magnetoe¬ 
lastic couplings B a /s and all phonon properties are known 
EH 33] * we can explicitly calculate the dissipation kernel 
Gff (t , t'). The induced magnetic field hi(t ) in Eq. Ill) is a 
linear functional of the phonon coordinates, see theSup- 
plemental Material m- Assuming that the phonons are 
in thermal equilibrium at temperature T, the hi(t) are 
Gaussian non-Markovian random processes. While the 
average field ( hi(t )) = hi(t) leads only to a small correc¬ 
tion to the external field H(t) (which we neglect), the 
fluctuating part Shift) = hi(t) — hi(t) is crucial for the 
thermalization and the emergence of a magnon conden¬ 
sate. As usual, the covariance of the induced magnetic 
field is related to the damping kernel via the fluctuation- 
dissipation theorem, (5h?(t)5hj (t')) = TG°)f (t,t'). 

The LLG 0 describes the precession of the spins in 
an effective magnetic field, which includes the effect of 
spin-spin interactions as well as thermal fluctuations and 
damping caused by the random magnetic field obeying 
the fluctuation-dissipation theorem. However, unlike the 
conventional stochastic LLG the damping is described by 
a tensor Gwhich is nonlocal in both space and 
time due to the nonlocal spin-spin interactions induced 


by the exchange of phonons and the memory effects of 
the thermal phonon bath. This non-Markovian damp¬ 
ing corresponds to a thermal noise field with a colored 
spectrum generated by the phonon dynamics. In con¬ 
trast to the white noise in the conventional LLG this col¬ 
ored noise can also be expected to accurately describe 
the physical situation on small time scales and when 
the phonon and spin dynamics have the same time scale 
[351 03] . Also note that the stochastic process with co- 
variance 0 is not ergodic due to the explicit dependence 
of Gon the spin variables Fortu¬ 

nately, in the experimentally used parallel pumping ge¬ 
ometry where the external magnetic field is of the form 
H(t) = [H 0 + Hi cos(2w p t)]e2 (where w p is the pumping 
frequency and the static part H 0 of the magnetic field 
is much larger than the amplitude Hi of the oscillating 
part) we can approximate Gy (t, t') by an ergodic process. 
In this geometry the ground state of the spin system is a 
saturated ferromagnet magnetized in z-direction. Writ¬ 
ing the spins as Sift) = S[e z +rrii(t)\ it is then reasonable 
to expect that |m,;(t) <C 1. Given the fact that for YIG 
S ss 14 is rather large, we may substitute SJt) —> Se z 
in Eq. ([2]). The damping kernel G ij{t,t') is then inde¬ 
pendent of the state of the spin system; i.e., it is ergodic. 
Another benefit of this approximation is that G ij(t,t') 
now depends only on the differences Jf?.; — Rj and t — t' 
so that it is convenient to work in momentum space and 
consider mfc(t) = )>A e~ lk ' R, mi{t). After this simplifica¬ 
tion, we can solve our LLG numerically without further 
approximation. Technical details of the numerial solu¬ 
tion can be found in the Supplemental Material |37| . In 
the rest of this Letter, we present our numerical results. 

To make contact with the experiments m we have 
chosen the parameters for our simulation such that they 
describe a typical experimental setup: we consider a thin 
stripe of YIG with thickness d, = 6.7 fi m at temperature 
T = 300 K in a bias magnetic field H 0 = 1710 Oe x p, 
and a pumping field Hi = 0.03 Hq, where /r = 2 hb is 
twice the Bohr magneton. The fields are assumed to 
be oriented parallel to the longest axis of the stripe, 
which we identify with the z-axis. Aligning the axaxis 
perpendicular to the stripe, the wave-vectors of the in¬ 
plane magnons are of the form k = k y e y + k z e z = 
\k\sin8ke y + |A;| cosine- and the corresponding long- 
wavelength magnon dispersion in the absence of the 
pumping field is [33,; [35] 

E k = [H 0 + p s k 2 + A(1 - f k ) sin 2 8 k \ 1/_ 

x [H 0 + p s k 2 + Af k ] 1/2 . (3) 

Here the energy A = 4npMs is determined by the satura¬ 
tion magnetization [33ll34ll45] 47 tMs = 1750 G, the spin 
stiffness is given by p s = JSa 2 = 5.17 X 10 -9 Oe cm 2 x n, 
and the form factor is ,33; iHO fk = [1 — e~^ d }/{\k\d). 
With the lattice spacing a = 12.4 A we obtain for the 
effective spin S = Msa 3 /p ~ 14.2. The magnon disper- 
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a) t = 15 ns b) t = 55 ns 



FIG. 1. (Color online). Snapshots of the stochastic time evo¬ 
lution of the magnon distribution rik(t)/n relative to the 
thermal distribution n^ 1 ', showing (a) the magnetoelastic hy¬ 
bridization, (b) the parametric instability, (c) the redistribu¬ 
tion of magnons in momentum space, and (d) the quasiequi¬ 
librium with the magnon condensate. 


sion assumes then its minimal value .Emin = 27r x 4.9 GHz 
for k = ±k m - ln e z with k m = 4.999 x 10 4 cm -1 . To fix 
the spin-phonon interaction and hence the damping ker¬ 
nel Gij(t,t') we need the effective ionic mass [5 3f55I[3T] 
M = 5.17g x a 3 /cm 3 , the velocities of the longitudi¬ 
nal and the transverse phonons, cm = 7.209 x 10 5 cm/s, 
c_l = 3.843 x 10 5 cm/s, and the magnetoelastic coupling 
tensor, which for the cubic lattice is of the form B a p = 
S a pB\\ + (1 — 5 a p)B±. For YIG, the relevant couplings 
are @3, 33,; B L ~ 2Ey = 6.96 x 10 6 erg x a 3 /cm 3 . For 
the pumping frequency we choose w p = 27T x 7.046 GHz, 
which is slightly above the ferromagnetic resonance at 
Efc = o = 27 t x 6.824 GHz. A convenient choice of the 
phonon polarization vectors e^\ appearing in Eq. |2]) can 
be found in Ref. [ 35 ] • All parameters in Eqs. 0 and 0 
are now fixed so that the results of our numerical simu¬ 
lations do not contain any adjustable parameters. 

In the following we present plots of the square of the 
transverse spin component nk(t) = \rrv^(t)\ 2 + |m|(t)| 2 . 
The average (n,k(t)) is directly proportional to the 
magnon number of spin wave theory [34 1 . In particular, 
in thermal equilibrium ( Uk{t)) = n j(, h oc T/Ek is pro¬ 
portional to the classical Rayleigh-Jeans distribution. In 
Fig. |T]we illustrate the main stages of the time evolution 
of rifc(t), assuming that the system is prepared in a fully 
polarized ferromagnetic state when the pumping field is 
switched on. The time evolution can be characterized 
by four distinct regimes: The population of the magnon 
spectrum begins on the energy surface where magnetic 
and elastic modes hybridize, as shown in Fig. |T] (a). 
Magnon-phonon hybridization in YIG has recently been 
discussed in Ref. [35]. Note that our non-Markovian de¬ 



FIG. 2. (Color online), (a) Logarithmic plot of the paramet¬ 
ric resonance. The blue solid line shows rik(t)/n for the 
parametrically unstable magnon with momentum k/k m i n = 
—2Ae y + 0.2 e z . For comparison, the red dotted line is the 
exponential growth predicted by linear spin wave theory as 
discussed in the main text. Note that in our nonlinear theory 
the growth saturates at some finite value close to the steady 
state predicted by S-theory (black dashed line). At longer 
times, the magnon distribution relaxes to a quasiequilibrium 
value, (b) Plot of the stochastic time evolution of the uniform 
mode \ml. =0 \ 2 /N 2 , which grows rapidly as soon as the growth 
of parametrically unstable modes saturates at t ss 50 ns, and 
saturates itself when the parametric magnons thermalize to 
quasiequilibrium around t ~ 100 ns. 


scription is crucial for obtaining the magnon-phonon hy¬ 
bridization [48]. At the next stage shown in Fig. 0(b) 
the growth of the magnon distribution is dominated by 
the parametric instability of magnons with energies Ek 
in the vicinity of the pumping frequency cu p . An example 
of a stochastic trajectory of the parametrically pumped 
magnon with Ek = w p is displayed in Fig. [2] (a). Note 
that linear spin wave theory predicts that in the regime 
of parametric instability the envelope of the magnon oc¬ 
cupation increases exponentially 133137], n k {t) oc e 2oikt 
with a\ = \V k Hi\ 2 — \Ek —w p | 2 . Here 14 is proportional 
to the so-called ellipticity of the spin waves and can be 
expressed in terms of the Fourier transform D’jf of the 
dipolar tensor as 14 = S(D^ X — E^ y )/(4E fe ). Because 
the growth coefficient ak depends on the ellipticity, the 
parametric pumping is most effective for elliptic magnons 
close to k z = 0. 

Because of interaction effects (which are nonpertur- 
batively taken into account in our approach) the expo¬ 
nential growth of the parametrically pumped magnons 
stops after approximately 40 ns. From Fig. [2] (b) we 
see that at this point the uniform longitudinal mode 
m| =0 starts to be populated, although the occupation 
always remains small. S-theory provides an estimate for 
the saturation value of the parametric magnons 03HI3] 
n k/ nt k ~ NSE k (\V k Hi\ - | E k - w p |)/(8TA), which is 
shown as a black dashed line in Fig. [2] (a) and agrees 
reasonably well with our simulation. In S-theory, this 
corresponds to a steady state, but our simulation reveals 
that at t « 50 ns the parametric magnons begin to relax 
to a quasiequilibrium value although the pumping field is 
still present. This is the third stage of the time evolution 
shown in Fig. [l] (c), where excited magnons generated 
via the parametric instability are scattered to the states 
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FIG. 3. (Color online), (a) Logarithmic plot of the stochas¬ 
tic time evolution of rik(t)/n \* of the condensing mode fc = 
fcminez (blue solid line), and of the parametrically unsta¬ 
ble mode with momentum fc/fc m i n = —2.4e y +0.2e z (red 
dashed line). Note that the growth of the condensate begins 
as soon as the parametrically excited magnons have relaxed 
to quasiequilibrium, (b) Time-averaged quasiequilibrium dis¬ 
tribution function (rife(t))/n*. h (blue solid line) on the axis 
k y = 0. The red dotted line shows the fit given in Eq. Q 
with [left = 0.995 -Emm, while the orange dashed line is the 
thermal part alone. The condensate contribution was mod¬ 
eled as a Gaussian centered at fc = fc m i n e 2 with standard 
deviations a z = 0.16 fc m i n and a y = 0.1fc m in, which are both 
smaller than the lattice spacing Afc = 0.2fc m j n . The inset 
shows a magnified view of the fit away from the condensate. 

with k y = 0. Eventually, in the fourth stage shown in 
Fig. |T] (d), the magnons accumulate at k z = fc m ; n to form 
the Rayleigh-Jeans condensate. The stochastic time evo¬ 
lution of the condensed mode is shown in Fig. [3] (a), with 
the primary parametric modes for comparison. The for¬ 
mation of the condensate is completed at t ss 200 ns; at 
later times, there are only small oscillations around a sta¬ 
ble quasiequilibrium state for all modes. This allows us to 
obtain thermal averages by averaging over time. Choos¬ 
ing the interval between 200 and 800 ns for the average 
we obtain the quasiequilibrium distribution 

(n.fc(t))/ri^. oc fi e ff) T Ti c 8(k T fcmin^z)- (4) 

The first part corresponds to a rescaled Rayleigh-Jeans 
distribution with effective chemical potential /i e s ss 




FIG. 4. (Color online), (a) Logarithmic plot of the decay 
of the condensate. The blue line denotes nfc(f)/nfc h at fc = 
fcmin e z , whereas the red dashed line is an exponential fit with 
the condensate lifetime r c ss 250 ns. (b) Plot of the stochastic 
time evolution of the uniform mode |m|_ 0 | 2 /./V 2 . As soon as 
the pumping is turned off at t = 100 ns the uniform mode 
starts to decay. 


F m i n , which is valid for all modes with the exception of 
the modes at the dispersion minima, for which one has 
to add delta-like contributions with weight n c and width 
comparable to the resolution Afc of the momentum grid. 
This proves that the classical spin dynamics indeed gives 
rise to Rayleigh-Jeans condensation of magnons at the 
minima of the energy dispersion. A plot of the time av¬ 
eraged distribution function as well as of the fitting func¬ 
tion 0 is displayed in Fig. [3] (b). Let us emphasize that 
the retardation of the dissipation kernel in Eq. 0 is not 
essential for the appearance of the condensate. However, 
if we use a phenomenological dissipation kernel with a 
shorter correlation time a condensate emerges only for 
unrealistic values of the magnon-phonon coupling and 
the pumping field. 

Finally, we show in Fig. [4] the time evolution of the 
system for the case that the pumping field is switched off 
at t = 100 ns. The condensate still forms in this regime, 
in agreement with the experiments mu (8j- However, 
at t « 250 ns it begins to decay, whereas the uniform 
mode decays as soon as the pumping field is turned off. 
Moreover, from the left panel it is clear that the decay 
of the magnon condensate is to a good approximation 
exponential, with characteristic decay time r c ~ 250 ns, 
in excellent agreement with the experiments mm. 

In summary, we have derived and solved a non- 
Markovian stochastic LLG which describes all stages 
of the nonequilibrium time evolution of the pumped 
magnon gas in YIG, including the magnon-phonon hy¬ 
bridization at early stages, followed by the parametric 
resonance regime and the emergence of a magnon con¬ 
densate, and the eventual decay of the condensate and 
the thermalization of the magnon distribution after the 
pumping has been switched off. Such a complete descrip¬ 
tion of all stages leading to the quasiequilibrium conden¬ 
sation of magnons in YIG has not been achieved before. 
Our results are in quantitative agreement with experi¬ 
ments and prove that the emergence of a condensate in 
the pumped magnon gas in YIG is a purely classical ki- 
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netic condensation phenomenon, similar to the conden¬ 
sation of classical light which has been observed in pho- 
torefractive crystals l30j . 

We thank A. A. Serga and D. A. Bozhko for discussions 
and the DFG for financial support via SFB/TRR 49. 
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SUPPLEMENTAL MATERIAL 

The following supplementary material contains techni¬ 
cal details of the calculations presented in the main text 
and some additional derivations. 

Stochastic Landau-Lifshitz-Gilbert equation from 
microscopic spin-phonon dynamics 


polarization A and velocity c\. The canonical momenta 
Pk\ and coordinates X^\ associated with the normal 
modes satisfy the commulation relations [X^xYk' v] = 
iN5k.-k’5\\i. With our normalization of the Fourier ex¬ 
pansion, the lattice displacement X(Ri) at lattice site 
Ri is given by 

X(Ri) = ^ E e lk RlX kxekx, (S6) 

k\ 


Here we derive the stochastic non-Markovian Landau- 
Lifshitz-Gilbert equation given in Eq. (1) of the main text 
starting from the quantum mechanical Hamiltonian for 
the coupled spin-phonon system in YIG. Our derivation 
follows the work by Rossi, Heinonen, and MacDonald 
m- The established Hamiltonian of the relevant spin- 
and lattice degrees of freedom in YIG consists of three 
parts [ S2] . 


H(t)=H m (t) + H e + H me , (SI) 


where the magnetic part T-L m {t) depends only on the spin 
degrees of freedom, T~L c {t) is the elastic (phonon) part, 
and 7~L me is the magnetoelastic coupling. Assuming that 
the system is exposed to a time-dependent external mag¬ 
netic field H(t), the magnetic part of the Hamiltonian 
is 

TMt) = ~\ 

ij a/3 

(S2) 

i 

where Si are quantum mechanical spin operators local¬ 
ized at sites Ri of a cubic lattice, the are the ferro¬ 
magnetic exchange couplings connecting nearest neigh¬ 
bors on a cubic lattice, and the matrix elements Dff of 

L J 

the dipolar tensor are (/jlb is the Bohr magneton), 


Dtf = (!-&„) 


( 2 Hb) 


2 r 


3M - 5 af} 


\Rij \ 3 L « « 


(S3) 


where R l: j = Ri — Rj and Rij = R,j/\R t j\. Defining the 
tensor with matrix elements K“' 3 = S a ^Jij + D“-' 3 , 
the spin Hamiltonian (S2) can alternatively be written 
as 


(t) = - \ E s ^ s i - E H (*) • 


(S4) 


The elastic part of our Hamiltonian is given by 


*.=sE 


N 


k X 


P-k\Pk\ M 2 v v 
2 M + Y kxX ~ kxXkX 


(S5) 


where N is the number of lattice sites and e^x are the po¬ 
larization vectors associated with the normal mode with 
momentum k and polarization A. Finally, the magnetoe¬ 
lastic coupling in YIG is to leading order in a gradient 
expansion given by (S2, S3 






(S7) 


a./3 


where the strain tensor is defined in terms of the com¬ 
ponents X a (r) of the phonon displacements in direction 

^•a: 

-.-a/3 = 1 [ 9X a (r) dXp{r) 

1 2 _ drp 8r a 

For YIG the magnetoelastic coupling tensor is of the form 

isiEa 



Pap $apP\\ A (1 ^a@)P -L- 


(S9) 


From the Hamiltonian (SI) it is straightforward to ob¬ 
tain the quantum mechanical equation of motion for the 
spin operators, 

Si(t) = Si(t) x [ff(t) + E X b^(0 

3 

+ * [Si{t) x Fi(t) - F t (t) x Si(t)] , (S10) 

where the vector operators F t (t) describe the effect of the 
magnetoelastic coupling on the spin-dynamics, 

W) = -4 E B a0 e a Xf(t)S?(t). (Sll) 


Oi(3 


The equation of motion (S10) should be complemented 
by the equations of motion for the normal modes of the 
phonons, 


Xkx = 


PkX 


M ’ 

Pkx = —M^kxXkx + MA k x(t ), 


where 


(512) 

(513) 


Ak x(t) = 


2 MS 2 


EE- 


-ik-Ri 


B, 


ol/3 


i ot(3 

xk af} -e- k xS?(t)S?{t). 


where M is the effective ionic mass of a unit cell and 
WfcA = c x \k\ is the dispersion of acoustic phonons with 


(S14) 
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Here fc aj g = k a ep + kpe a . Taking an additional time 
derivative of Eq. (S12) we obtain 


(dt + ^kxjXkX — Ak\(t), 
which has the general solution 

X k xft) = X kx (0) cos(w fcA t) + P , k r X ^ sin(wfc A t) 

1 V 1 U>k\ 


(S15) 


/o 




= -A fc A(0) cos(wfcAt) + sin(wfcA^) 

Mu>k A 

+ ^UaW _ f 4 ^, cos[u;fcAft ~ *' )] dA k \(t') 


u ; 


kX 


LU 


kX 


dt' 


(S16) 


In the last line, we have integrated by parts and intro¬ 
duced the notation 


A fc A(0) = X feA (0)-^^. 


UJ 


(S17) 


k A 


At this point we replace the spin operators by classical 
vectors, ignoring the fact that different spin components 
do not commute. The magnetoelastic field F t {t) can then 
be written as the following functional of the spins, 

F i (t) = h i {t) + 6h i {t)- [ dt'VG lj (t,OS i (t'), (S18) 
J o 


where the part h t (t) of the induced magnetic field is in¬ 
dependent of the phonon dynamics, 

Ht) = ^iEE ^B a pB^ Y, {t)Sj{t)Sj ft) 

ap p,u j 


"nX 

kX 


,ik-(Ri~Rj 


i {J^otfi ' ^kx){^ , iiv ‘ &—kx) 


M “L 


(S19) 


while the part Shift) depends on the initial conditions of 
the phonons coordinates and momenta at time t = 0, 

Shift) = --±-Y B a pe a Sf(t) Y e ik R 'k a0 ■ e k x 


cx.fi 


kX 


A’fcA(O) cos(wfcAi) + -rr^~ sm{w k \t) 
Muj k \ 


Finally, the damping kernel is given by 

Y B ^ B p- s i^) s W) 


(S20) 


x ~^Y elk(R ' Rj) ( k w ‘ e k\)(kp„ ■ e_ k \) 

k\ 

COs[w fc A(f - t’)] 


M “lx 


(S21) 


Let us now assume that the (shifted) initial values X k \ (0) 
and Pk\(0) of the normal modes are classical Gaussian 
random variables with zero average and variance given 
by the classical equipartition theorem with temperature 
T. The averages then vanish, (X k \(0)) = 0 = (PkxfO)), 
while the second moments are 

(X fc A(0)AW(0)> = JV<S fci _ fc ,<S aa'vt^, (S22a) 

M “tx 

(Pkx(0)Pk'X'm=NS k ,-k'S X yMT, (S22b) 

<*hA(0)JW(0)>=0, (S22c) 

where (...) denotes the Gaussian averaging over the clas¬ 
sical thermal phonon distribution. Note that this implies 


(Shift)) = 0, (S23) 


and that the covariance of the random fields Shfft) is 
related to the matrix elements of the damping kernel 
(t, f 7 ) via the fluctuation-dissipation theorem, 

(6hf(t)6h^t'))=TGff(t,t'). (S24) 


With hi ft) = hi ft) + Shift), we finally obtain the 
stochastic Landau-Lifshitz-Gilbert equation with non- 
Markovian damping given in Eq. (1) of the main text, 


Ht) = S^t) x 


- Sift) X 


H(t) + hi(t) + Y^ij s j(t) 

j 

f dt' Y, t')Sj(t'). (S25) 

J o j 


Assuming that the phonon dynamics is much faster than 
the magnon dynamics, and that the damping kernel is 
local in the spatial indices and diagonal in the spin labels, 
we may approximate 


Gffft, t') » 2 'ySijS^Sft - t’). (S26) 

In this approximation Eq. ( |S25[ ) reduces to the usual form 
of the Markovian LLG 


Sift) = Sift) x 


H{t) + h i {t) + Y^i3Sj(t)-'rHt) ■ 

j 

(S27) 


However, because the characteristic phonon energies in 
YIG have the same order of magnitude as the magnon 
energies, the approximation (S26) is not sufficient for a 
realistic description of the interplay between magnon and 
phonon dynamics in YIG. 


Numerical solution of the LLG 

Reduction of the non-Markovian random process to 
Gaussian white noise 


In order to numerically solve the non-Markovian 
stochastic Landau-Lifshitz-Gilbert equation (S25) it is 
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convenient to reduce it to a Markovian equation contain¬ 
ing only Gaussian white noise by introducing additional 
variables. The only approximation necessary is to ne¬ 
glect the dependence of the induced field and the damp¬ 
ing kernel on the state of the spin system by replacing 


Si(t) Se z in their respective definitions, Eqs. (S20) 
and (S21). Then it is convenient to define the Fourier 


transforms 

8h k (t) = Y, e~ lk R 'Sh t (t) = Y th kX (f), (S28) 

i X 

G k (t, 0) = Y e~ lk Rij Gij(t, 0) = Y cos (w k \t), 

i A 

(S29) 

where 

8h k \ (f) — — ^ ' B (XZ s OL k c 


e k \ 


A fcA (0) cos (w k \t) + PkX ^ sin (ui k \t) 

M u k a 

(S30) 

ir d d (k az ■ e k \) {kp z • e_ k \) 

G kX = B az B 0z - Mw 2 a52 -■ ( S31 ) 

The corresponding momentum space version of the 


fluctuation-dissipation theorem (S24) is 


{5hl x (t)5hl, x ,(t')) = N5 k - k '5\yTG al3 


kX 


x cos \u> k \(t — t')] . (S32) 

In the following we will neglect the average field hi(t) 
since it only gives rise to a small correction of the external 
field, which we take as 


H(t) = [H 0 + Hi cos(2w p f)]< 


(S33) 


Also, we will use that for thin YIG films the dipolar ten¬ 
sor is diagonal in the spin components [S4]. Expanding 


the spins in the LLG (|S25j) as 

Si{t) = S 


1 

N y 6 

k 


£■ 


ik-Ri 


m k (t) 


(S34) 


yields the equation of motion for the Fourier modes 

m k (t), 

rh k = [L fe + P (t)] m k + e z xY j^fcA(i) 


-SG k x / dt' cos [uj kx {t - t')\ rn k (t') 




q x SK q m q + Y [dh qX (t) 


—SG q a / dt 1 cos [uj q \(t — t')\ rn q {t') 

J 0 


(S35) 


Here we introduced the matrices 


Li, = 


0 

-wS 



and 


0 1 0> 

P(f) = Hi cos(2w p t) (—100 
0 0 0 ; 


(S36) 


(S37) 


Lfc describes the linear oscillations of the spins in terms 
of the frequencies 

^ = Ho + S(Kf =0 -Kn, (S38) 

^ = Ho + S(Kf =0 -K™), (S39) 

whereas P(f) contains the parametric pumping. Note 
that the nonvanishing eigenvalue of the matrix h k gives 
the magnon dispersion of linear spin wave theory, 


Ek = y/u%u>l (S40) 

The simple harmonic time dependence of the damping 


tensor (S29) now allows us to convert this non-Markovian 


equation of motion into a system of Markovian ones |S51 
S(i . To this end we introduce the new set of variables 

U k A =SG k A / dt 1 cos [u} kx {t - t')\ rn k (t') 

Jo 

— SG k \m k . (S41) 

They satisfy the second order equations of motion 


U kX — U) kx U k A SG kX TTl k 
with initial conditions 

u k a(0) = -SG k \m k (0), u k A (0) = 0. 


(S42) 


(S43) 


With the definition ( S41| ) of the new variables the spin 
equation of motion (S35) becomes 


rn k = [L fe + P (t)} m k 

+e z x Y [ Sh k\(t) ~ u kx - SG kX m k 


-Y 

AT 


m k _ 


N ^ ---9 
Q 


SK q m q 


~ ^ A(^-) HqX SG q \Tn q 

X 


(S44) 


The remaining task is to simulate a Gaussian ran¬ 
dom process with covariance given by the fluctuation- 
dissipation theorem (S321. As will be shown below, such 


a process can be modeled with a driven oscillator equa¬ 
tion of motion [S5tiS7l 


(dt + 2gwfe A i9t + Wfe A ) 8h k a = 2^jgTu kX b kX f kX {t), 

(S45) 
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where the driving term fk\(t) is a complex Gaussian 
white noise process with vanishing average and covari¬ 
ance 


= N6 k k'8xx'8(t ~ t')- 


(S46) 


The coupling vector is given by 


bkX ^ ^ ^q-^o 


ik n 


e k x 


xfMuJkxS 


(S47) 


so that b kx b®_ kx = G kx - We also added a damping term 
with strength g < 1 which is necessary to obtain a finite 
covariance function and can be interpreted as the damp¬ 
ing of the phonon subsystem. In the end we have to take 
the limit g —» 0 to reproduce the Fluctuation-Dissipation 


theorem (S321 from the driven oscillator (S45). To stay 
consistent, a corresponding term —2gujkx‘u>kx should also 
be added on the right-hand side of the equation of motion 


(S42) for Ukx- 


We have now succeded in our task to get rid of the non- 
Markovian terms and the colored noise, at the expense of 
having introduced two new sets of equations of motion. 
However, it turns out that we can combine these two 
sets into a single one since they are both linear oscillator 


equations and the spin equation of motion (S44) depends 
only on the difference 


Vkx = Sh k x - u k x , 


(S48) 


which also obeys a driven oscillator equation. Therefore 
we are left with the two coupled sets of stochastic differ¬ 
ential equations 


m k = [L k + P(i)] m k 


+e z x 


£[■ 


^ k\ 


N 


rn k _ q x l SK q m, 


Y[ 


VqX — SG q \m q 
VkX SGkX^'k 


(S49a) 


Vkx = -ZguikxVkx ~ ukx 

+2 jgTi4 x b k xfkx(t). (S49b) 

with complex white noise fkx(t ), and the constraint 

Vkx( 0) = S<Gkxm k ( 0), u fcA (0) = 0 (S50) 

on the initial conditions of the auxiliary variables. 


Numerical implementation 

A plot of the magnon and phonon dispersions for YIG 
with parameters mentioned in the main text is given in 



k / k min 


FIG. SI. (Color online). Dispersions of the magnon and 
phonon modes of a thin YIG stripe with thickness d = 6.7 /xm 
in an external magnetic field Ho = 1710 Oe x /x, for k parallel 
and perpendicular to the in-plane magnetic held. 


Fig. |S1| Especially note that both phonon branches 
cross the magnon dispersion relatively close to the min¬ 
imum of the magnon dispersion. Therefore phonon and 
magnon dynamics occur at the same time scale and we 
cannot approximate the damping tensor (S21) by white 
noise as in the conventional Markovian LLG. However we 
will neglect the longitudinal phonon branch entirely and 
only consider the two transverse branches. This can be 
justified with the fact that the damping tensor (S21) is 
proportional to the square of the inverse sound velocity. 
Therefore the faster longitudinal phonons are consider¬ 
ably less important for the magnon relaxation than the 
transverse phonons. A plot of the corresponding energy 
surface E k = uik± is shown in Fig. [S2j together with 
the contour satisfying the parametric resonance condition 
Ek — w p - Note however that the parametric resonance is 
only effective for elliptic spin waves which have momenta 
close to k z = 0 in thin YIG films. The phonon damping 
coefficient necessary for the noise simulation will be set 
to g = 10 -3 . 

For the explicit numerical solution we use the stochas¬ 
tic Heun scheme [EE1 ES] to solve the system (S49) of 
stochastic differential equations on an equidistant mo¬ 
mentum grid with increment Afc = 0.2/c m j n . As we are 
mainly interested in the long-wavelength regime we work 
with a truncated momentum space with 125 x 189 = 
23625 points, such that \k y \ < 12.4fc m i n and \k z \ < 
18.8 fc m ; n . The integrations are discretized as 





l 


£- 

(S51) 


where we defined the effective lattice size N e g = 
47r 2 /(aAfc) 2 ss 2.6 x 10 7 . The corresponding delta dis¬ 
tributions are replaced according to NSkk' —• y N e g5kk' 
since A NSkk' = 1 has to hold. As all momentum 
integrations are convolutions, they are performed using 
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ky / bnill 


FIG. S2. (Color online). Momentum space contours of the 
energy surfaces of parametric resonance Ek = to p (blue solid 
line) and transverse magnetoelastic hybridization Ek = tafcx 
(red dashed line). 


a fast Fourier transform algorithm jSlOj . Gaussian ran¬ 
dom numbers are created via the polar Box-Muller algo¬ 
rithm ISllj from random numbers distributed uniformly 
between 0 and 1. We calculate the stochastic time evolu¬ 
tion of mfc(t) from t = 0 up to t = 800 ns in increments 
of At = 10 -4 ns, starting from a fully polarized ferro¬ 
magnetic state with 0 ) = 0 for all momenta k. 


Noise driven harmonic oscillator 


This appendix will be devoted to simulating complex 
colored noise satisfying 

<*(*)) = 0, (S52) 

(x(t + r)x*(t)} = K cos(f It), (S53) 

via an underlying white noise process. To this end we 
consider an (under-)damped, complex harmonic oscilla¬ 
tor driven by Gaussian white noise, i.e., 

(d? + 2'yd t + n 2 )x = Af(t), (S54) 


with 


</(*)> = (S55) 

if{t + T)f*(t)) =5(t). (S56) 


The general solution of Eq. 
situation is given by 


(S54|) in the underdamped 


(S57) 


with the Green function 


G(t) = Q(t)e 


-7 1 


sin ^ yji I 2 — 7 2 tj 

VW-'I 2 


(S58) 


and where Xh(t) denotes the solution of the homogeneous 
equation. As we aim to simulate a stochastic process 
with vanishing average, we set Xh(t) = 0 by choosing 
x(0) = 0 = x(0). In frequency domain, we then obtain 


c(w) = / dte lut x(t) = AG(w)/(w) (S59) 


where 


G M = 


1 


(S60) 


W-oj 2 - 2iju] 

For the random force, we have 

(f(t + T)f*(t))=6(r) 

. [ dwi [ ^2 -iun(t+ T )+iu 2 t 

J 27T J 2tt 

X (/( Wl )r( W2 )), (S6i) 

and therefore 

(/(wi)/*(w 2 )) = 2ttS(oji - w 2 ). (S62) 

With the above, we can calculate the correlation function 


(x{t + T)x*(t)) = \A\ 2 J 
= l^| 2 / 


dw 

27T f 

dw 


-|GM| 2 


(S63a) 


(ft 2 


■ uj 2 Y + 4y 2 w 2 


(S63b) 


l^l 2 „-7ld 


4 7 ft 2 ~ cos(V^ 2 -7 2t ) 

■ ^J_^ sin( v /0 2 -7 2 M) 


(S63c) 

For our application, we are mainly interested in the limit 
of vanishing damping. To this end we define the dimen¬ 
sionless damping constant g = 7 /ST, so that 


{x{t + r)x*(t)) = -J 4 Le _fln|r| 


4 g n 3 


+ 


g n 


V 1 - a 2 M 


cos ^\/l — < 7 2 ftr^ 
sin (\/l - 3 2 |nr|) 


(S64) 


Next, we demand that 


lim(ir(t + t)x*U)) = K cos(ftr), (S65) 
g ->0 


yielding 


\A\ = s/4gWK. 


(S 66 ) 


x(t) = x h (t) + A dt'G(t - 
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Parametric instability 


Thermal equilibrium 


In this appendix we will derive the parametric instabil¬ 
ity from linear spin wave theory. Therefore we consider 
only the first line of the equation of motion (S49a) for 
the spin deviation, which explicitly yields 


TOfc = \(jj v k + H X cos(2u;pf)] m k , (S67a) 

fn v k = — [uj k + Hi cos(2w p f)] m k . (S67b) 


The parametric instability is usually derived for the com¬ 
plex normal mode of the system without pumping, given 
by 


<t>k = 





(S68) 


In terms of this normal mode the two linearized equation 
of motions (S67) become the single complex equation 


fih ^ \E k T UfoHi cos(2w p t)] (j) k 

— 2iV k H\ cos(2ui p t)(j)*_ k , (S69) 


where U k = y/l +4|I4| 2 and V k = (u v k - wj£)/(4 E k ) = 
S(D k x — D k v )/(4E k ) is proportional to the ellipticity of 
the spin wave. Next let us move to a reference frame 
rotating with the pumping frequency by setting 


= e lulpt tp k (t). 


(S70) 


Inserting this ansatz into the equation of motion (S69) 
yields 


^ p k — ^ T U k H\ cos(2w p t)] Vk 

-iV k Hi (e 4iWpt + 1) <f*_ k . (S71) 

To study the parametric instability it is now sufficient 
to only retain the resonant term and drop all oscillating 
contributions in the above equation [S3]. This leaves us 
with 


Vk = ~i{E k - u> p )vk - iV k Hnp*_ k . (S72) 

Taking an additional time derivative turns this into 

Vk = -(\E k -u p \ 2 - \V k Hi\ 2 )v k , (S73) 

showing that the system becomes unstable for 

\E k -u p \<\V k Hi\. (S74) 

This is the parametric instability [S3] . Note that it can 
only occur for elliptic spin waves, which are character¬ 
ized by nonvanishing V k . It is also readily seen that the 
unstable modes diverge as Vk oc e at , with 


(S75) 


To describe thermal equilibrium in a ferromagnetic 
state we may safely assume that the x and y compo¬ 
nents of the spin deviation vector m.j = (S', — Se z )/S 
are small compared to unity. Then we may expand 

m* = Jl- (mf) 2 - (mf) 2 - 1 « [(ml) 2 + (mf) 2 ] . 

(S76) 

Retaining only terms quadratic in mf and ml in the 
Hamiltonian (S4) and dropping the pumping field leads 
to 




—SNH 0 — E ] 


Hj 




= -SNH 0 - 


n 

H 0 + SJ2 

n 

s 2 


TZZ 

^in 


ij oi=x,y 


-SKg“ to“to“ 


zz 

fc =0 


^E E 

k a=x,y 

(S77) 

As we are now restricting ourselves to a noninteracting 
system, it is again convenient to introduce the normal 
modes (S681. In terms of them the above Hamiltonian is 
given by 

H m = -SNHo - ^Kf =0 + ^ E (S78) 


which is just a sum of ground state energy and an en¬ 
semble of noninteracting, classical harmonic oscillators. 
In thermal equilibrium each of them has to satisfy the 
equipartition theorem 

2 T 

<■<&</>*)= N5 kk . —, (S79) 

and therefore the thermal equilibrium value of the square 
of the transverse spin component n k = |m£| 2 + |ro^| 2 is 
given by 

27 1 

n? = (nk) = Ny/l + 4|Vfc| 2 ——(S80) 
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